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Simulation  of  inviscid  compressible  multi-phase 
flow  with  condensation 


Philip  KeUeners  f 


Condensation  of  vapours  in  rapid  expansions  of  compressible  gases  is  investigated.  In 
the  case  of  high  temperature  gradients  the  condensation  will  start  at  conditions  well  away 
from  thermodynamic  equilibrium  of  the  fluid.  In  those  cases  homogeneous  condensation 
is  dominant  over  heterogeneous  condensation.  The  present  work  is  concerned  with  de- 
velopment of  a simulation  tool  for  computation  of  high  speed  compressible  flows  with 
homogeneous  condensation.  The  resulting  flow  solver  should  preferably  be  accurate  and 
robust  to  be  used  for  simulation  of  industrial  flows  in  general  geometries. 


1.  Introduction 

A substance  below  its  critical  temperature  can  be  present  in  either  gaseous  or  liquid 
phase,  depending  on  the  pressure,  and  is  referred  to  as  a vapour.  Vapours  present  in 
a mixture  of  gases  and  vapours,  when  subjected  to  expansion  can  condensate  and  form 
liquid  droplets.  This  phenomenon  is  observed  in  e.g.  aircraft  tip  vortices  and  in  industrial 
flows  like  steam  turbines.  Condensation  in  flows  of  gas  mixtures  at  high  speed  has  been 
investigated  by,  among  others;  Wegener  (1969),  Hill  (1966),  Campbell  (1989),  Schnerr  et 
al.  Schnerr  (1996),  Dohrmann  (1989),  Mundinger  (1994),  Adam  (1996)  and  van  Dongen 
et  al.  Luijten  (1999),  Luijten  (1998),  Prast  (1997)  and  Lamanna  (2000).  Expansion  in 
nozzles  of  gases  to  supersonic  speeds  has  often  been  used  to  investigate  the  physics  of 
condensation.  Condensation  in  the  flow  around  airfoil  sections  and  in  steam  turbines  has 
been  investigated  to  a large  extent.  At  the  University  of  Twente  a numerical  tool  has 
been  developed  to  simulate  transonic  flows  with  condensation  in  confined  geometries. 
The  solver  operates  on  the  basis  of  a finite  volume  method  using  unstructured  meshes. 
It  has  been  observed  that  results  obtained  with  the  solver  are  very  sensitive  to  accurate 
shock  prediction,  and  fine  shock  resolution  in  the  flow  field,  especially  in  cases  of  strong 
interaction  between  the  gasdynamic  shock  and  the  condensation  process.  The  focus  of  the 
present  work  is  to  improve  the  accuracy  and  robustness  of  the  flow  solver  by  improving 
solid  wall  boundary  treatment  and  spatial  reconstruction  for  simulations  with  second 
order  spatial  accuracy. 


2.  Physics  of  Condensation  during  Rapid  Expansion 

Below  its  critical  temperature,  a fluid  can  be  in  gaseous  or  liquid  phase.  The  thermo- 
dynamical region  of  coexistence  of  vapour  and  liquid  in  equilibrium  in  bulk  substances 
is  given  by  the  Clausius-Clapeyron  relation.  For  rapid  expansions  of  vapour  this  coexis- 
tence region  is  passed  without  the  fluid  attaining  equilibrium.  The  vapour  is  saturated 
expressed  by  the  saturation  ratio; 
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Figure  1.  Expansion  of  air-water  mixture  in  nozzle-S2,  partial  vapour  pressures  and 

saturation  ratio  s (dashed) 


where  n„  and  n„,eg  are  the  number  of  moles  of  the  vapour  in  the  actual  mixture  and 
the  mixture  in  equilibrium,  respectively.  For  pressures  as  high  as  atmospheric  pressure 
equation  2.1  can  be  replaced  by  the  more  conventional  expression: 


Pv,eq 

where  is  the  actual  partial  vapour  pressure  in  the  mixture  and  Pv,eq  is  the  equilibrium 
partial  vapour  pressure  at  the  same  thermodynamic  conditions.  s,p„  and  Pv,eq  are  plotted 
in  figure  1 for  the  expansion  of  an  air-water  mixture,  as  functions  of  temperature. 

The  curve  for  the  coexistence  of  liquid  and  vapour  {Pv,eq)  divides  the  plane  given  by 
pressure  and  temperature  into  two  regions.  For  pressures  higher  than  the  coexistence 
pressure  the  substance  under  consideration  (water)  is  present  in  liquid  form  while  in 
equilibrium  state.  For  lower  pressures  the  substance  will  be  present  as  water-vapour 
while  in  equilibrium.  The  expansion  in  the  nozzle,  as  given  by  the  partial  vapour  pressure 
Pv,  is  so  quick  that  the  water-vapour  will  not  immediately  condense  under  equilibrium 
conditions.  This  is  the  case  as  the  characteristic  time  of  the  gasdynamic  flow  is  much 
smaller  than  the  time  needed  to  form  the  first  onsets  of  the  new  liquid  phase.  For  the 
nozzle  flow  at  hand  this  is  due  to  the  high-cooling  rate  in  the  nozzle.  So  the  water-vapour 
expands  further,  driving  the  air-water-vapour  mixture  well  away  from  equilibrium,  as 
indicated  by  the  saturation  ratio  which  attains  values  as  high  as  40.  In  case  s > 1 the  fluid 
is  said  to  be  super-saturated.  Formation  of  small  liquid  clusters,  nuclei,  at  high  super- 
saturation is  the  first  stage  of  the  condensation  process  that  starts  in  order  to  reestablish 
equilibrium.  On  these  newly  formed  nuclei,  the  super-saturated  vapour  condenses  as  a 
second  step  until  equilibrium  is  reached,  the  process  of  droplet  growth.  This  can  be  seen 
in  figure  1 by  the  decrease  in  partial  vapour  pressure  and  the  increase  in  temperature. 
With  so  much  liquid  already  present,  the  remainder  of  the  condensation  process  due 
to  droplet  growth  takes  place  very  close  to  thermodynamic  equilibrium  of  the  water, 
indicated  by  the  saturation  close  to  one,  for  the  remainder  of  the  expansion. 
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The  process  described  above  is  the  process  of  homogeneous  condensation.  In  contrast, 
the  process  of  heterogeneous  condensation  is  droplet  growth  on  foreign,  already  present 
particles.  However,  for  high  cooling  speeds  and  consequentially  high  super-saturation 
the  number  of  new  nuclei  formed  exceeds  any  realistic  number  of  foreign  particles  by 
several  orders  of  magnitude.  For  the  present  applications  only  homogeneous  condensation 
processes  are  of  importance. 

An  effect  resulting  from  condensation  is  the  release  of  latent  heat  during  the  con- 
densation process,  indicated  by  the  increase  in  temperature  during  condensation  in  the 
previous  example  of  nozzle  flow.  The  latent  heat  L is  defined  as 

L = K-hi  (2.3) 

with  hv  and  hi  the  enthalpy  of  the  gaseous  and  liquid  phase,  respectively.  The  latent 
heat  is  the  enthalpy  needed  to  evaporate  a unit  mass  of  liquid.  It  is  a material  property. 


3.  Nucleation  and  Droplet  Growth  models 

From  the  previous  treatment  of  condensation  during  rapid  expansion,  it  is  observed 
that  the  condensation  process  consists  of  two  consecutive  stages.  The  first  one  is  the 
formation  of  liquid  nuclei,  nucleation,  the  second  one  is  the  condensation  of  vapour 
molecules  on  the  already  present  nuclei,  making  these  grow  in  the  process  of  droplet 
growth.  In  the  following,  a brief  treatment  is  given  of  the  physical  background  of  both 
nucleation  and  droplet  growth  together  with  the  presentation  of  the  models  used  to 
simulate  these  processes. 

3.1.  Nucleation 

Under  super-saturation  vapour  molecules  can  condense  on  liquid  already  present  or  on 
foreign  objects.  In  the  absence  of  both,  the  vapour  molecules  can  form  small  clusters.  As 
a result  of  the  clustering,  an  additional  phase  needs  to  be  formed,  the  interface  between 
the  liquid  inside  the  cluster  and  the  gas  outside  of  the  cluster.  The  interface  can  be  re- 
garded as  infinitely  thin.  At  the  interface  there  is  surface  tension.  Thus  the  creation  of  an 
interface  requires  energy.  For  very  small  clusters,  the  surface  effects  are  dominant  over 
volume  related  effects.  As  a result  the  formation  of  the  interface  represents  an  energy 
barrier  in  the  formation  process  of  the  nucleus.  If  the  energy  involved  in  the  clustering 
of  the  vapour  molecules  is  less  than  required  in  the  formation  of  the  interface  surface, 
the  cluster  will  disintegrate  immediately  following  its  formation.  Therefore  at  near  equi- 
librium conditions,  super-saturations  close  to  one,  it  is  highly  unlikely,  that  a realistic 
number  of  stable  nuclei  in  a volume  at  macro  scale  will  be  formed  although  clusters 
are  constantly  formed  and  falling  apart.  The  previous  notion  results  in  formulation  of 
formation  enthalpy  of  critically  sized  stable  nuclei,  and  the  number  densities  in  which 
these  stable  nuclei  are  likely  to  come  into  existence.  The  creation  enthalpy  of  one  liquid 
droplet  under  equilibrium  conditions  is  given  by  Dohrmann  (1989): 

AG  = 4irr^(T  — nfcTln(s)  (3.1) 

where,  AG  is  the  change  in  Gibbs  enthalpy,  r is  the  radius  of  the  droplet,  cr  is  the  droplet 
surface  tension  in  a plane  with  no  curvature,  n is  the  number  of  molecules  in  the  droplet, 
k is  the  Boltzmann  constant,  and  T is  the  temperature.  In  accordance  with  classical 
nucleation  theory,  we  note  that  a stable  nucleus  will  be  created  when  the  function  for 
AG  attains  a local  maximum.  AG  at  this  maximum  is; 
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AG  = |irr*  V (3.2) 

with  r*,  the  radius  of  this  critical  sized  droplet  being  given  by: 

^ piKTlnis) 

where  pi  is  the  density  of  the  liquid  in  the  droplet,  and  is  the  gas  constant  of  the 
vapour.  The  number  of  droplets  created  is  given  by  the  nucleation  rate.  For  the  nucleation 
rate  several  models  are  available,  Classical  Nucleation  Theory  CNT-model  (Wegener 
1969)  and  the  Internally  Consistent  Classical  Theory,  ICCT-model  (Luijten  1998).  In 
the  present  work,  the  CNT-model  is  applied: 

/ 16rr 

Y 3 mp‘}RlT^\n'^{s) 
where  is  the  density  of  the  vapour  and  m is  the  mass  of  one  vapour  molecule. 

3.2.  Droplet  Growth 

The  droplet  growth  model  to  be  applied  depends  on  the  regime  in  which  droplet  growth 
takes  place.  For  pressures  of  the  order  of  atmospheric  pressure  droplet  growth  is  based 
on  a balance  between  condensation  of  vapour  molecules  onto  droplets,  and  evaporation 
of  vapour  molecules  from  the  droplet.  For  pressures  1 to  2 orders  of  magnitude  higher, 
droplet  growth  is  diffusion-controlled.  In  the  present  work,  high  pressure  effects  are  not 
taken  into  account.  Therefore  the  Hertz-Knudsen  droplet  growth  model  (Hill  1966)  can 
be  applied.  The  droplet  growth  rate  is: 


dr  _ a / Pl£__\  /■OK') 

dt  Pi  \ ^ ’ 

where  a is  an  accommodation  coefficient  usually  taken  equal  to  one,  Td  is  the  droplet 
temperature  in  the  present  work  equal  to  the  surrounding  gas  temperature  and  Ps,r  is 
the  partial  super-saturated  vapour  pressure  over  a curved  radius  of  curvature  r given  by: 


Ps,r  — Vv,eq  6Xp 


2a  \ 

piRuTrhi) 


(3.6) 


where  rt^i  is  the  Hill  droplet  radius.  This  radius  is  computed  as  an  averaged  radius  over 
the  complete  population  of  droplets  at  hand.  This  radius,  of  course,  depends  on  the 
distribution  function  of  the  droplets. 

4.  Description  of  the  Liquid  Phase 

Until  now  only,  the  creation  and  growth  rate  of  a single  droplet,  and  the  creation  rate, 
nucleation  rate,  of  new  droplets  has  been  treated.  Next  to  that,  the  size  and  form  of 
individual  droplets  and  the  size  and  distribution  of  the  droplet  population  need  to  be 
computed.  As  the  number  of  droplets  in  the  flow  cases  at  hand  typically  range  between 
10^^  and  10^^  number  of  droplets  per  unit  mass,  it  is  impossible  to  describe  individual 
droplets.  One  alternative  approach  is  to  define  classes  of  droplets  of  similar  radius,  the 
class-model.  This  approach  naturally  extends  to  a description  of  the  distribution  of  the 


Compressible  multi-phase  flow  with  condensation 


53 


Figure  2.  Liquid  droplets  in  control  volume 

droplets.  For  the  nozzle  flow  problems  at  hand,  the  typical  range  of  radii  of  the  droplets 
will  be  between  2.10“^°  and  1.10“'^.  This  would  require  an  impractical  high  number  of 
droplet-classes,  making  the  class-model  computationally  very  expensive.  Hill’s  method 
of  moments  is  computationally  less  intensive  at  the  cost  of  some  of  the  resolution  of  the 
droplet  distribution.  As  for  the  problems  at  hand,  the  total  amount  of  liquid  generated 
and  the  strong  interaction  between  the  gasdynamic  flow  and  the  condensation  process  is 
of  primary  interest,  the  loss  of  resolution  of  the  droplet  distribution  is  no  severe  penalty. 

4.1.  Hill’s  Method  of  Moments 

Consider  a control  volume  as  depicted  in  figure  2.  In  the  control  volume  a mixture  of  an 
ambient  gas,  and  a vapour  in  its  gaseous  and  its  liquid  state  is  present.  A first  assumption 
is  that  all  liquid  is  present  in  the  form  of  spherical  droplets.  Conservation  of  the  mass  in 
the  control  volume  gives: 


}VI 


M ^Ma  + M^  + Mi  (4.1) 

where  M is  the  total  mass  in  the  control  volume,  Ma  is  the  mass  of  the  ambient  gas, 
Mv  is  the  mass  of  the  condensible  substance  in  gaseous  form  and  Mi  is  the  mass  of  the 
same  substance  in  liquid  form.  The  ambient  air  is  assumed  to  be  permanently  above  its 
critical  temperature.  The  dimensionless  liquid  mass  firaction  g is  defined  as; 


Ml 

M ~ 


pV 


Vi 


(4.2) 


With  the  final  notation  emphasis  is  placed  on  the  next  assumption;  the  liquid  mass 
present  in  the  control  volume  is  regarded  to  be  incompressible.  Therefore  pi  is  taken 
constant.  The  liquid  mass  fraction  is  a function  of  time.  The  density  of  the  complete 
mixture;  gas,  vapour  and  liquid  p,  is  dependent  of  time,  as  are  the  volume  occupied  by 
the  liquid  V/  and  the  total  control  volume  V.  This  gives; 


9{t)=Pi- 


Vi{t) 


p{t)v{ty 

The  volume  occupied  by  the  liquid  Vj  is  given  by  the  following  integral: 


Vi{t)  = j -■nr^{t,r)J{T)V{T)dT  (4.3) 

0 

where  J(t)  is  the  nucleation  rate  per  unit  volume  and  V{t)  is  the  size  of  the  control 
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volume,  at  the  moment  of  nucleation.  The  volume  occupied  by  the  droplets  is  given 
by  ^TTr^{t,r).  The  function  r{t,T)  describes  the  radius  of  the  droplet,  from  its  time  of 
creation,  r(ro,  tq)  where  this  radius  is  equal  to  the  critical  radius  as  given  in  equation  3.3, 
r{To,To)  — r*,  up  to  the  present  radius  of  the  droplet  at  time  t.  Under  the  integral  the 
product  J(to)U(to)  acts  as  a weight  function  for  the  contribution  of  the  liquid  volume 
of  the  droplets  created  at  time  tq  to  the  total  amount  of  liquid  formed  at  time  t.  The 
liquid  mass  fraction  now  becomes: 

J ^nr^(t,T)J(T)V{T)dT 

0 


If  the  considered  lump  of  mass,  liquid  and  gas-vapour  mixture,  is  at  rest,  or  flows  along 
a stream-line,  the  fraction  is  fixed.  This  is  identical  to  the  assumption  that  the 

droplets  do  not  move  relative  to  the  surrounding  gas  mixtiu-e.  This  is  known  as  the 
no-slip  condition  in  Hill’s  Method  of  Moments.  This  can  be  expressed  as; 


p{t)V{t) 


dM  = 0=^M  = p{t)V{t)  = p(t)V{t) 

The  liquid  mass  fraction  can  be  rewritten  (Hagmeijer  2001); 

g{t)=Pi  J ^7rr^(t,T)^^dT. 

0 

However,  the  function  r(t,  r)  is  not  readily  available  in  closed  form.  The  droplet  growth 
rate  ^ is  known  in  closed  form.  Therefore  g{t)  is  differentiated  with  respect  to  time: 


dg  d 

It  ~ 


P(t) 


dr 


-/I 


7r3r^(t,r) 


dr{t,r)  J(t) 
dt  p{t) 


, 4 ,J(t) 

dr  + p,-nr 


— n ^ *^(^)  I o. 


dr{t,r)  J(t) 

; T-rdr. 

dt  p{t) 


In  this  analysis,  it  is  assumed,  that  at  the  initial  time  t — 0 there  is  no  liquid  present 
in  the  control  volume.  A third  assumption  is  that  the  present  droplet  growth  rate  is 
independent  of  the  present  radius  of  the  droplet.  This  assumption  is  certainly  not  true 
for  very  small  droplets  because  of  the  dependency  of  the  surface  tension  on  the  droplet 
radius,  however  the  assumption  is  valid  for  droplets  with  larger  radii. 


%^%irit,r))  (4.4) 

This  implies  that,  the  present  value  of  the  droplet  growth  rate,  is  independent  of  its 
history.  So  it  can  be  written: 


dr 

dt 


(4.5) 
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In  this  case,  it  is  allowed  to  take  the  variation  of  the  radius  with  time  out  of  the  integral: 


J 1 4/3nrHt,r)^dr. 

0 

Careful  inspection  of  the  right-hand  side  of  the  relation  above  shows,  that  the  original 
integral  has  reappeared;  but  with  the  power  of  the  radius  reduced  by  one,  and  multiplied 
by  the  power  of  the  function  r and  the  droplet  growth  rate.  Following  this  observation, 
it  is  helpful  to  employ  the  following  definition  of  a liquid  moment: 

t 

Qn{t)  = J r"(t,r)^^dT. 


Qs  is  closely  related  to  the  liquid  mass  fraction: 


/..N  ^ V 4 dQsit) 


The  variation  of  this  liquid  moment  Qn(t)  with  time  is  very  much  similar  to  expression 


dQn 

dt 


' p{t)  dt 


t 

{t)  j r"-i(t,r) 


JJj) 

P{r) 


dr 


which  is  a recurrent  relation  for  the  differentiation  of  Qn-  First  consider  the  formulation 
of  the  relevant  liquid  moments: 


Inspection  of  the  last  moment,  Qo,  shows  that  there  is  no  longer  a functional  dependence 
of  t in  the  integral.  So  the  differentiation  of  this  last  liquid  moment  can  be  expected  to 
be  different  from  the  higher  liquid  moments; 

dt  dt  yj  p{t)  ^ j p[t)' 
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This  last  result  is  fortunate,  as  it  shows,  that  the  logical  expansion  defined  by  the  recur- 
rence relation  for  the  liquid  moments,  is  not  continued  for  n < 0.  So  there  is  no  closure 
problem  for  the  solution  of  the  set  of  moments.  The  system  of  ordinary  differential  equa- 
tions, to  determine  the  liquid  moment  Qs  now  becomes: 


dQs 

dt 

dQi 

dt 

dQi 

dt 

dQo 

dt 


...dr 

dr 

, J(t)  dr  , 

Pit) 


In  this  system  the  variables  represent  the  following  quantities; 

• Qa  is  proportional  to  the  sum  of  all  droplet  volumes, 

• Q2  is  proportional  to  the  sum  of  all  droplet  surface  areas, 

• Qi  is  proportional  to  the  sum  of  all  droplet  radii, 

• Qo  is  the  present  number  of  droplets. 

Now  it  is  assumed  that  the  droplet  growth  rate  ^(t)  is  given  only  by  the  present  local 
thermodynamic  and  chemical  state,  and  not  by  the  spatial  gradients  of  r in  the  flow 
domain.  In  this  case,  = 0.  By  addition  of  the  product  of  vector  Q and  the  continuity 
equation  for  mass,  the  system  of  ordinary  differential  equations  above  can  be  rewritten 
in  strong-conservation  form: 


dpQs  dpQsUj 

dt  dxi 

dpQ2  dpQ2Uj 

dt  dxi 

dpQi  dpQiUj 

dt  dxi 

dpQo  dpQpUj 

dt  dxi 


= r*V(t)-p3j(t)pQ2 

= r-^J(t)  + 2^{t)pQi 

fiv 

= r* J(t)  -h  ^^{i)pQo 
= J{t) 


6 p 

Where  Xi  and  Ui  are  the  spatial  coordinates,  and  the  velocity  components  in  spatial 
direction  respectively,  r*,J  and  ^(t)  are  given  by  the  laws  formulated  for  the  creation 
and  growth  of  a single  droplet.  The  above  set  of  equations  gives  the  integral  properties 
of  the  droplet  distribution.  So  details  of  individual  droplets  are  not  available.  To  be  able 
to  compute  a radius  needed  for  droplet  growth  laws,  Hill  defined  (Hill  1966)  an  averaged 
radius  as: 


(4.6) 
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Reason  for  application  of  Q2  in  calculation  of  r^i,  is  the  dominance  of  surface  effects  in 
droplet  growth,  hence  the  application  of  the  integral  surface  quantity  Q2- 


5.  Thermodynamic  State  of  Gas-Liquid  Mixtures 

The  change  of  phase  from  vapour  to  liquid  releases  latent  heat  to  the  surrounding 
mixture.  Due  to  the  condensation,  the  fractions  of  the  gases  in  the  mixtures  change. 
Both  processes  result  in  a change  of  the  thermodynamic  state  of  the  mixture.  To  model 
these  changes  regard  the  enthalpy  of  the  mixture: 


Mh  — Maha  + hdyh'u  + Mihi 

where  h = H/M  is  the  specific  enthalpy.  By  application  of  the  previously  defined  dimen- 
sionless mass  fraction  5 = 

h — (1  5max)^a  T (f/max  Q')^v  "h  Qhl- 

Application  of  the  definition  of  the  latent  heat  L = hy— hi  and  the  definition  Cp  = )p  ^ 
results  in: 


— (1  5max)Cp„  -t-  {gmax  fl)^v  d"  Oippy  Qrp  )•  (^-1) 

Under  the  assumption  ^ ~ ^ ^ similar  expression  can  be  derived  for  the  isochoric 

specific  heat  coefficient  Cy  = {§f)yg- 


^ — (f  5max)Cua  d"  (,9max  9)^v  d"  9{^v  )•  (®'2) 

With  these  relations  for  the  specific  heat  coefficients,  the  gas  constant  Rmix  aud  the 
dimensionless  ratio  'fmix  can  be  computed  similar  to  the  case  of  an  ideal  gas,  Rmix  (9,L)  = 
Cp  — Cy  and  7mix(5,  L)  = Both  R and  7 are  now  functions  depending  on  the  mass 
fraction  g and  the  relation  for  the  latent  heat  L.  Derived  quantities  like  pressure  p and 
speed  of  sound  c can  be  shown  (Mundinger  1994)  to  be: 

p = pR(g,l)T  c^  = 'i{g,L)?-^. 


6.  Model  of  Inviscid  Flow  with  Condensation 

Based  on  Reynolds  numbers  for  typical  flow  problem  of  interest  the  flow  is  assumed 
to  be  inviscid.  The  Euler  equations  are  used  to  describe  the  flow.  Together  with  Hill’s 
Method  of  Moments,  the  models  for  nucleation  and  droplet  growth  and  an  equation  of 
state  they  form  a closed  set  of  equations: 

^ J (h'dv+  j E{(p')  ds  = j W{(^')dv  (6.1) 

V S=9V  V 

where  <p'  is  the  vector  with  the  conserved  quantities,  F{4i')  the  flux  vector,  and  W {<!>')  is 
the  source  term: 
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where  p,u,p  are  the  density,  the  velocity  and  the  pressure  of  the  mixture.  E = e + 
and  H = E-\-^  are  the  total  energy  and  the  total  enthalpy  of  the  mixture.  Qi  is  the  i-th 
moment  in  Hill’s  Method  of  Moments.  The  thermal  equation  of  state  reads: 

p = pR(g,L)T  (6.3) 

and  the  caloric  equation  of  state: 

e = c^{g,L)T.  (6.4) 


7.  Finite  Volume  Discretization 

With  the  definition  of  the  control  volume  averaged  4>i: 

' Vi 

equation  6.1is  semi-discretized  as: 

j=i 

n 

The  summation  of  discrete  fluxes  is  computed  using  an  edge-based  data 

i=i 

structure  as  described  in  Jameson  (1986),  Barth  (1989)  and  Barth  (1994).  For  every 
edge  the  indices  of  the  2 control  volumes  connected  by  the  edge  are  stored,  as  well  as 
the  three  spatial  components  of  the  surface  normal  vector  of  the  surface  between  the  two 
control  volumes.  The  length  of  the  surface  normal  vector  is  equal  to  the  magnitude  of  the 
surface  area.  This  edge-based  approach  allows  for  a cell-centered  or  a vertex-centered 
use  of  the  original  mesh,  without  significant  changes  to  the  flux  computation  algorithm. 
The  fluxes  are  computed,  by  solving  a local  one-dimensional  Hiemann  problem  for  every 
surface  using  an  approximate  Riemann  solver.  This  is  implemented  in  the  flow  solver  for 
edges  oriented  in  every  possible  direction  in  three-dimensional  space.  Two-dimensional 
problems  with  meshes  in  two  independent  directions  can  be  regarded  as  a subclass  of  the 
full  three-dimensional  problem.  The  same  holds  for  quasi  one-dimensional  or  truly  one- 
dimensional flow  problems.  With  proper  definition  of  additional  boundary  surfaces  in  the 
case  of  quasi  one-dimensional  flows,  all  flows  can  be  computed  using  the  same  flow-solver 
regardless  of  the  mesh  being  fully  three-dimensional,  two-dimensional  or  quasi  or  truly 
one-dimensional.  The  preprocessor,  to  the  flow  solver,  generating  the  edge-based  data 
structure  from  the  meshes  needs  only  to  insert  zeros  for  the  edge-related  surface  normals 
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in  the  directions  not  relevant  for  the  particular  mesh.  To  date,  pre-processors  have  been 
written  for  meshes  consisting  of  line  elements  for  one-dimensional  meshes,  triangular  ele- 
ments for  unstructured  two-dimensional  meshes,  quadrilaterals  for  structured  monoblock 
two-dimensional  meshes,  and  tetrahedral  elements  for  three-dimensional  meshes.  How- 
ever, the  edge-based  data  structure  allows  simple  extension  to  multiblock  meshes,  meshes 
constructed  by  hexahedra  and  hybrid  meshes.  This  requires  only  the  pre-processor  for 
the  particular  element  (s)  to  be  developed,  and  minor  additions  to  the  domain  boundary 
treatment  in  the  flow  solver.  The  flow  solver  has  successfully  completed  simulations  for 
one-,  two-  and  three-dimensional  flows.  The  great  advantage  of  this  approach  is  that 
almost  all  testing  during  development  can  be  done  for  flow  problems  in  one  and  two- 
dimensions  decreasing  development  time,  as  formulation  of  one  and  two-dimensional  flow 
problems  and  generation  of  one-  and  two-dimensional  meshes  is  much  less  time  consum- 
ing than  in  the  full  three-dimensional  case.  There  is  a penalty  in  the  case  of  computation 
of  one-  and  two-dimensional  flows.  For  the  one-dimensional  case  the  momentum  fluxes 
in  the  two  spatial  directions  perpendicular  to  the  main  spatial  direction  are  computed 
but  not  used.  In  the  case  of  the  Euler-equations  this  would  result  in  66%  increased  com- 
putational expense.  In  the  two-dimensional  case  this  increased  computational  expense 
is  25%.  In  the  case  of  the  computation  of  the  Euler-equations  and  Hill’s  momentum 
equations,  this  computational  overhead  drops  to  28%  in  the  one-dimensional  case,  and 
12%  in  the  two  dimensional  case.  However,  the  number  of  grid  points  for  flow  prob- 
lems in  two-dimensions  and  one-dimensions  is  one  and  two  orders  lower  respectively.  So 
the  computational  overhead  is  completely  negligible  next  to  the  enormous  advantage  of 
simple  and  faster  testing  in  lower  dimensions. 

Eigenvalue-analysis  Kelleners  (2001)  of  the  complete  system  of  Euler-equations  with 
the  Hill  Momentum  equations  has  shown  that  the  additional  eigenvalues  due  to  the  Hill 
equations  are  all  real  and  have  value  u.  So  the  liquid  moments  Qo  -Qs  are  convected 
downstream  along  streamlines,  and  no  new  acoustic  waves  are  introduced.  The  computa- 
tion of  the  flux  is  identical  to  the  case  of  flow  without  condensation,  with  only  additional 
transport  equations  for  the  liquid  moments. 

The  presence  of  the  source  terms  in  the  Hill  momentum  equations  results  in  the  com- 
plete system  of  partial  differential  equations  being  a stiff  system.  Wishing  to  use  an  ex- 
plicit time-stepping  method  as  a relaxation  method  to  compute  a steady  state,  the  fact 
that  the  system  is  stiff  would  lead  to  an  unacceptable  small  time-step  requirement  due 
to  the  source  terms.  To  circumvent  this  time-step  restriction  a fractional  time-stepping 
method  is  used  as  proposed  by  Oran  md  Boris,  Oran  (1987)  the  differential  equation  of 
the  form: 


+ fW  = w[4>) 


the  solution  procediure  is  split  into; 


d<t>* 

dt\ 

dti 


-f{4>n 

w{4>*) 


(7.1) 

(7.2) 


where  the  time-step  operator  for  the  first  step  can  be  any  conventional  explicit  operator, 
e.g.  Euler-forward  or  Runge-Kutta.  The  time  step-operator  for  the  second  time  step 
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needs  to  be  able  to  cope  with  the  possible  large  magnitude  of  the  source  w{4>*),  (see 
Oran  1987;  Mundinger  1994) 

The  two  following  subsections  highlight  two  topics  which  require  detailed  attention  in 
the  case  of  development  of  a spatial  second-order  accurate  finite-volume  method.  For 
other  topics  like,  flux  splitting  or  flux  limiting,  see  Jameson  (1995)  and  Liou  (1996). 

7.1.  Higher  Order  Spatial  Reconstruction 

To  achieve  good  accuracy  on  meshes  with  moderate  vertex  densities  the  computational 
method  must  be  at  least  second  order.  Whereas  computation  of  first-order  spatial  ac- 
curacy allows  for  the  storage  location  of  the  cell-averaged  value  to  be  located  anywhere 
inside  the  control  volume,  this  point  should  preferably  be  located  at  the  geometric  center 
of  gravity  for  higher-order  spatial  accuracy.  Use  of  the  center  of  gravity  circumvents  se- 
vere penalties  to  the  applicable  CFL-number  of  individual  control  volumes.  Higher-order 
reconstructions  of  cell-averaged  data  should  conserve  the  cell-averaged  value.  In  case  of 
second-order  reconstruction  it  can  be  proven  that  use  of  the  center  of  gravity  satisfies 
this  requirement: 


u(x)  =Ucg  + -^  (x-  x^g)  (7.3) 

where  Ugg  is  the  cell-averaged  value  positioned  at  the  center  of  gravity.  u{x)  integrated 
over  the  control  volume  gives; 


^ J u{x)dv  = y J Uggdv+^  {x-Xgg)dv  (7.4) 

V V K “ 

The  first  term  on  the  right-hand  side  gives  the  cell-averaged  value.  The  second  term  on 
the  right  hand  side  is  zero.  The  gradient  can  be  moved  in  front  of  the  integral.  The 
remaining  integral  is  the  static  moment  about  the  center  of  gravity,  and  necessarily  zero. 
For  vertex-centered  control  volumes,  the  movement  of  the  cell  storage  location,  from 
cell-vertex  to  center  of  gravity  of  the  control  volume,  is  largest  for  control  volumes  at  the 
physical  boundaries  of  the  domain. 

7.2.  Linear  Reconstruction  at  Domain  Boundaries 
The  assumption  of  linear  reconstruction  in  the  control  volume  in  the  case  of  second- 
order  spatial  accuracy,  requires  special  attention  at  the  domain  boundaries.  At  solid  wall 
boundaries  the  pressure  part  of  the  flux  needs  to  be  calculated.  Question  is,  how  the 
pressure  must  be  integrated  along  the  solid  wall  boundary  surface  for  the  resulting  flux 
integral  along  the  complete  boundary  of  the  control  volume  to  be  of  second-order  spatial 
accuracy.  It  is  assumed  that  fluxes  vary  linearly  over  the  control  volume.  Reconstruction 
of  the  flux  at  the  solid  wall  boundary  is  therefore  similar  to  linear  reconstruction  from 
the  cell  averaged  state  of  the  conservative  variables  throughout  the  control  volume.  The 
latter  requires  the  computation  of  gradient  the  ^ at  the  center  of  the  control  volume. 
Green-Gauss  reconstruction  is  often  used  for  computation  of  this  gradient  Barth  (1994); 
Aftosmis  (1995).  Green-Gauss  reconstruction  is  easily  implemented  using  the  edge-based 
data  structure.  Starting  point  is  Gauss’  divergence  theorem: 

j ^i<i>)dv  = j 4>ds 
V s 


(7.5) 
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Figure  3.  Median  dual  mesh  control  volume  sharing  sides  with  a physical  boundary 
which  is  rewritten  in  discrete  form: 

i 

As  an  example  the  linear  distribution  of  <j>  along  a domain  boundary  for  a control  volume 
being  the  median  dual  to  a two-dimensional  triangular  element  is  given,  see  figure  3. 
The  volume  V is  the  surface  area  associated  with  vertex  0 in  the  case  of  the  median  dual 
mesh  is: 

^ = ^llaxb|| 

To  evaluate  equation  7.6,  the  outward  pointing  surface  vectors  S1...S4,  and  the  mean 
values  of  the  quantity  need  to  be  determined. 

1 

s_i  = -a  X cz 

Ah  1 ^ 

s_2  = (gg  - gS)  X e_z 
s_3  = (gO  - -a)  X e_z 
s^  — 2“  ^ ^ 

where  = e^,  x gj,.  4>i  is  computed  similar  as  would  be  the  case  for  a median  dual  mesh 
volume  embedded  entirely  in  the  physical  domain.  The  value  for  (j>i  at  surfaces  «2  and  S3 
is  interpolated: 


0Sg  = *^2) 

At  the  surfaces  s^  and  S4  it  is  assumed  that  there  exists  a linear  distribution  of  (j)  along 
the  surface  of  the  following  form: 


<j)  - uj4>o  + (1  — oj)4>i 
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The  allowable  value  of  w is  between  zero  and  one.  Substitution  of  and  4>i  in  equation 
7.6  ultimately  gives: 


X 


((3ai  — |)«^o  +(§  ~ 3a>)</>i 


-4>2) 


1 

Q-x 


((—3a;  + f )<^o 


+^i 


+(— I + 3o;)<?I>2) 


1 

O-x  ^y~^y^x 


To  compute  the  value  of  u>,  the  gradient  V(<^)  is  also  computed  assuming  a linear  distri- 
bution of  (j>  in  the  triangle  given  by  the  vectors  a and  b. 


(j){x)  - <j)o  = a{x  - xo)  -h  j3{y  - yo) 


The  gradient  for  this  linear  distribution  is: 

vw=(i)=(;) 

Using  the  data  in  vertices  1 and  2,  4>i,^2  ^ system  can  be  derived,  which  can  be  solved 
for  ^ ^ ^ . The  solution  of  this  system  is: 

f a \ _ f ay  \~^  f 4>i  - 4>o  \ 

\ 0 J \ b:,  by  J \ (f>2- 00  ) 


Both  expressions  for  ^(0)  are  equal  to  one-another.  Formulating  this  identity,  and  noting 
that  it  must  be  true  for  any  value  of  0o,0i<02,  it  can  be  shown  that: 


w = 


5 

6 


Returning  to  the  assumption  that  the  wall  pressure  flux  is  linearly  distributed  adong  the 
domain  boundary  surfaces,  in  a manner  similar  as  any  reconstructed  scalar  0,  the  weights 
I and  I can  be  applied  to  compute  the  integral  of  the  pressure  over  the  domain  boundary 
surfaces.  E.g.  for  surface  §4  this  becomes: 


1 , 

PS,S4  = (gP0+  gP2)s4- 


It  should  be  stressed  that  the  value  of  u depends  on  the  type  of  element  at  the  domain 
boundary,  (e.g.  triangle,  tetrahedron  or  brick)  and  the  manner  in  which  the  control 
volume  is  defined  (e.g.  cell-centered  or  vertex-centered  median  dual  mesh). 


8.  Resiilts 

The  simulation  tool  described  above  has  been  used  to  compute  solutions  to  flow  prob- 
lems of  internal  and  external  flows,  adiabatic  flows  and  flows  with  condensation.  The 
cases  presented  below  have  been  computed  using  the  following  layout  of  the  solution 
algorithm.  Node-centered  finite-volume  scheme,  using  median  dual  mesh  cells.  Explicit 
time  stepping  of  the  homogeneous  equations  as  described  in  equation  7.2  using  the  second 
order  Heun-method.  For  second  order  spatial  resolution,  application  of  a least-squares 
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algorithm  for  computation  of  spatial  gradients  in  the  flow  field.  Reconstruction  of  limited 
second-order  states  at  the  control  volume  interfaces  using  minimum-maximum  restric- 
tions formulated  by  Barth  (1989),  and  a limiter  function  by  Venkatakrishnan  as  presented 
in  Aftosmis  (1995).  Fluxes  computed  using  flux  splitting  according  to  Eberle.  Compu- 
tation of  the  time-updates  due  to  the  source  terms  as  mentioned  in  equation  7.2  using 
a fractional  time  stepping  method  as  presented  in  Mundinger  (1994).  Inflow  boundaries 
conserving  entropy  and  total  enthalpy  at  infinity  upstream  for  a general  equation  of  state, 
outflow  boundaries  maintaining  downstream  static  pressure  in  case  of  subsonic  outflow, 
or  simple  first-order  extrapolation  of  flow  quantities  in  case  of  supersonic  flow.  At  solid 
walls  application  of  linear  pressure  distribution  along  the  outward  facing  surfaces,  as 
presented  in  section  7.2.  Heat  effects  due  to  condensation,  are  taken  into  account  by  in- 
clusion of  both  the  gaseous  and  liquid  phase  into  the  energy  conservation  equation.  This 
results  in  strong  two-way  coupling  of  the  condensation  process  and  the  gasdynamics. 
All  computations  related  to  the  equation  of  state,  or  material  properties  have  been  pro- 
grammed into  separate  library  routines.  These  libraries  are  called  upon  by  the  routines 
solving  the  conservation  equations.  This  produces  additional  computational  overhead, 
but  should  allow  for  quick  implementation  of  different  equations  of  state.  In  the  present 
form  of  the  flow  solver  only  an  equation  of  state  for  an  ideal  gas  is  implemented. 

8.1.  Condensation  in  Nozzle  Ba-120 

Flows  with  condensation  in  nozzle  Ba-120,  designed  by  Bartlma,  were  investigated  both 
experimentally  and  numerically  by  Schnerr  and  co-workers  (see  Mundinger  1994).  Flows 
for  two  different  humidities  are  presented  in  figure  4.  The  stagnation  conditions  for  both 
flows  in  figure  4 are: 

so[%]  To[K]  po  [pa] 

42.1  298.8  100900 
49.3  297.1  100900 

The  expansion  of  the  air-water-vapour  mixture  gives  rise  to  condensation  in  the  flow 
field  downstream  of  the  nozzle  throat.  This  is  seen  by  the  steep  rise  in  the  nucleation 
rate.  Following  the  creation  of  the  nuclei  is  the  process  of  droplet  growth  indicated  by 
the  rise  of  the  liquid  mass  fraction  —2—.  The  release  of  latent  heat  to  the  flow,  induces 

* Smax 

a shock  in  both  cases.  However  the  shock-strength  in  the  case  of  the  higher  humidity  is 
much  larger,  and  as  a result  the  nucleation  pulse  is  terminated  abruptly.  In  both  cases, 
the  continuing  expansion  due  to  nozzle  divergence  results  in  condensation  of  almost  all  of 
the  water-vapour  at  nozzle  station  x = .15[m].  In  the  Mach-isoline  plot  of  the  flow  with 
higher  humidity,  the  reflection  of  the  curved  shock-wave  from  the  nozzle  wall  is  nicely 
visible. 


8.2.  Condensation  in  Vortex  Flow 

To  study  the  process  of  condensation  in  vortical  flow,  a very  slender  delta  wing  has  been 
placed  inside  a tube,  see  figure  5.  At  the  inlet  of  the  tube  a mixture  of  air-water-vapour 
flows  in  at  Mach  1.55.  The  humidity  of  the  mixture  at  stagnation  condition  is  as  low  as 
1.1%.  The  wing  induces  a vortex  in  the  flow  field.  In  the  vortex  core  there  is  considerable 
loss  of  total  pressure,  clearly  visible  in  figure  6.  The  vortex  spirals  downstream  in  the 
tube  in  a helical  form.  In  the  low  pressure  region  in  the  vortex  on  the  upper  side  of  the 
wing,  super-saturation  occurs,  resulting  in  large  nucleation  rate,  see  figure  7.  Note  that 
the  region  with  the  highest  nucleation  rates  remains  restricted  to  the  front  half  of  the 
wing  close  to  the  apex.  This  is  a consequence  of  the  droplet  growth  lowering  the  levels  of 
super-saturation  considerably  once  the  first  liquid  droplets  are  created  in  great  numbers. 
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Figure  4.  Flow  with  condensation  in  nozzle  Bal20,  flow  from  left  to  right.  Humidity  at  stag- 
nation condition  for  left  figure:  so  = 42.1%,  right  figure  so  = 49.3%.  Mach  iso-lines  in  upper 
figures,  only  iso-lines  for  M > 1 are  drawn.  Nucleation  rate  J,  dimensionless  pressure  and 
liquid  mass  fraction  — 2_  in  lower  figures.  Dashed  lines  indicate  center  line,  and  position  of 

* 9max  ' *■ 

nozzle  throat 

The  dimensionless  liquid  mass  fraction  — ^ is  plotted  in  figure  8.  On  the  second  half  of 
the  upper  side  of  the  wing,  almost  all  of  the  water  is  present  in  liquid  form.  Beyond  the 
trailing  edge  of  the  wing,  some  evaporation  of  the  liquid  water  occurs,  but  the  main  part 
of  the  water  is  convected  downstream  in  liquid  form.  As  the  vortex  is  diffusing  over  the 
cross-sections  of  the  tube,  so  is  the  cloud  of  water  droplets,  as  can  be  noted  from  the 
similar  shapes  of  the  iso-lines  of  both  the  total  pressure  in  figure  6 and  the  liquid  mass 
fraction  in  figure  8 


9.  Conclusions 

A brief  treatment  has  been  given  on  condensation  in  high  speed  flows.  Because  of  large 
temperature  gradients,  the  condensation  can  start  well  away  from  thermodynamic  equi- 
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Figure  5.  Layout  of  the  delta  wing,  which 
produces  the  vortex.  For  clarity  only  one 
half  of  the  tube  wall  is  shown.  Flow  from 
upper  right  corner  to  lower  left  corner. 


Figure  6.  Total  pressure,  pt 


Figure  7.  Nucleation  rate,  J plotted  on  a „ o t • -j  r 

log-scale  Figure  8.  Liquid  mass  fraction,  g. 


librium,  leading  to  fog-like  clouds  consisting  of  many  small  droplets.  Basic  nucleation 
and  droplet-growth  models  have  been  presented.  A thorough  derivation  of  the  integral 
description  of  the  liquid  phase,  Hill’s  Method  of  Moments  has  been  given.  The  com- 
plete model  describing  inviscid  compressible  multi-phase  flow  with  condensation  and  its 
finite-volume  discretization  have  been  given.  Two  aspects  of  the  implementation  of  the 
finite-volume  method  for  achieving  second-order  spatial  accuracy  have  been  highlighted. 
Results  presented  for  a two-dimensional  flow  problem,  and  a fully  three-dimensional  flow 
problem  have  been  computed  using  one  and  the  same  flow  solver.  The  benefit  of  develop- 
ing only  one  single  flow  solver  is  clearly  felt.  From  the  case  of  higher  humidity  in  nozzle 
Ba-120,  a strong  interaction  between  gasdynamics  and  condensation  can  already  be  seen. 
Higher  humidities  will  lead  to  unsteady  flow.  Both  experimental  and  numerical  data  is 
already  available  from  Schnerr  and  co-workers.  A future  extension  of  the  flow  solver  to 
take  into  account  unsteady  flow  should  be  made.  For  flows  with  condensation  in  vortices 
simulations  should  be  performed  using  the  full  Navier-Stokes  equations,  to  validate  the 
assumption  of  inviscid  flow,  made  in  this  work.  Improvements  must  be  made  with  respect 
to  accuracy  and  robustness  of  the  solution  algorithm,  to  make  results  of  the  simulation 
tool  less  dependent  on  grid  quality,  but  the  last  problem  is  universal  to  CFD-methods. 
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